In [2]:
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from libpysal.weights import Queen
import esda
from esda.moran import Moran, Moran_Local
import contextily as ctx
from splot.esda import moran_scatterplot, plot_moran, lisa_cluster,plot_moran_simulation
import libpysal as lps
import matplotlib.pyplot as plt
import plotly.express as px
from shapely import wkt
from shapely.wkt import loads
import branca.colormap as cm
import folium
import mapclassify

df = pd.read_csv("merged_MAR16.csv")

gdf = gpd.GeoDataFrame(df, geometry=df['wkt'].apply(loads))
gdf.set_crs(epsg=3857, inplace=True)

gdf = gdf.dropna()

invalid_geometries = gdf[~gdf.is_valid]
print(f"Found {len(invalid_geometries)} invalid geometries")
gdf['geometry'] = gdf.geometry.buffer(0)

gdf['snap_rate'] = gdf['snap_rate'].astype(float)
gdf['chd_pct'] = gdf['chd_pct'].astype(float)
gdf['lowaccess_pct'] = gdf['lowaccess_pct'].astype(float)
gdf['lowaccess_li_pct'] = gdf['lowaccess_li_pct'].astype(float)
gdf['snap_rate'] = gdf['snap_rate'].astype(float)
gdf['total_pop'] = gdf['total_pop'].astype(int)

gdf['centroid'] = gdf.geometry.centroid
centroids_geo = gdf['centroid'].to_crs(epsg=4326)


gdf['longitude'] = centroids_geo.x
gdf['latitude'] = centroids_geo.y

gdf = gdf[["geometry", "snap_rate", "chd_pct", "lowaccess_pct", "lowaccess_li_pct", "geoid20", "SPA_NAME", "pop_ages_65_older_pct", "latino_pct", "unemployed_pct"]]

Found 186 invalid geometries
In [3]:
spa_colors = {
    'West': 'red',
    'San Fernando': 'blue',
    'South': 'green',
    'San Gabriel': 'purple',
    'Antelope Valley': 'orange',
    'South Bay': 'yellow',
    'Metro': 'brown',
    'East': 'gray'
}

index_colors = {
    0: '#ffffcc', # light yellow
    1: '#ffeda0', # light orange
    2: '#fed976', # slightly darker orange
    3: '#feb24c', # orange
    4: '#fd8d3c', # dark orange
    5: '#fc4e2a', # reddish orange
    6: '#e31a1c', # red
    7: '#bd0026', # dark red
    8: '#800026', # burgundy
    9: '#400026', # dark burgundy
    10: '#000000' # black
}

spa_colors

gdf = gpd.GeoDataFrame(gdf, geometry=gdf["geometry"])
gdf.set_crs(epsg=3857, inplace=True)
Out[3]:
geometry snap_rate chd_pct lowaccess_pct lowaccess_li_pct geoid20 SPA_NAME pop_ages_65_older_pct latino_pct unemployed_pct
0 POLYGON ((-118.298 34.263, -118.297 34.263, -1... 0.000000 5.192139 0.024308 0.007821 101110 San Fernando 20.926756 26.955656 7.790493
1 POLYGON ((-118.277 34.260, -118.278 34.260, -1... 2.936857 4.704301 0.465066 0.037391 101122 San Fernando 21.829971 8.693563 9.590236
2 POLYGON ((-118.278 34.256, -118.278 34.253, -1... 2.701243 5.896306 0.000000 0.000000 101220 San Fernando 16.345877 42.028152 12.653688
3 POLYGON ((-118.287 34.256, -118.287 34.253, -1... 2.894356 5.398263 0.000000 0.000000 101221 San Fernando 19.515442 32.055378 4.291146
4 POLYGON ((-118.286 34.256, -118.286 34.255, -1... 6.915629 5.410280 0.000000 0.000000 101222 San Fernando 16.346153 46.118233 11.773151
... ... ... ... ... ... ... ... ... ... ...
2480 POLYGON ((-118.510 34.187, -118.508 34.187, -1... 44.843048 4.046243 0.000000 0.000000 980024 San Fernando 16.788321 15.328467 0.000000
2481 POLYGON ((-118.248 33.847, -118.247 33.847, -1... 6.180470 4.300235 0.503704 0.152189 980025 South Bay 26.277372 4.866180 13.592233
2484 POLYGON ((-118.281 33.767, -118.281 33.768, -1... 0.000000 4.518828 0.000000 0.000000 980031 Unknown 11.006289 34.591194 0.000000
2487 POLYGON ((-117.972 34.037, -117.972 34.037, -1... 21.279613 4.543835 0.032253 0.014417 980035 San Gabriel 10.917722 77.215187 10.263929
2488 POLYGON ((-118.079 34.043, -118.079 34.043, -1... 0.000000 5.430255 0.128709 0.046114 980036 Unknown 13.530928 47.680412 5.956113

2467 rows × 10 columns

In [4]:
classifier = mapclassify.NaturalBreaks.make(k=5)

# Age score
gdf['age_pct_score'] = gdf[['pop_ages_65_older_pct']].apply(classifier)
gdf[['pop_ages_65_older_pct', 'age_pct_score']].head()

# Hispanic score
gdf['hisp_pct_score'] = gdf[['latino_pct']].apply(classifier)
gdf[['latino_pct', 'hisp_pct_score']].head()

# Unemployment score
gdf['emp_pct_score'] = gdf[['unemployed_pct']].apply(classifier)
gdf[['unemployed_pct', 'emp_pct_score']].head()

# Food Index: Age + Hispanic + Unemployment scores 
gdf['priority_index'] = gdf['age_pct_score'] + gdf['hisp_pct_score'] + gdf['emp_pct_score']

gdf['priority_index'].unique()
Out[4]:
array([ 6,  5,  7,  4,  3,  8,  2,  1,  0,  9, 10])
In [7]:
chd_pct_min = gdf['chd_pct'].min()
chd_pct_max = gdf['chd_pct'].max()

snap_rate_min = gdf['snap_rate'].min()
snap_rate_max = gdf['snap_rate'].max()

lowaccess_li_pct_min = 0
lowaccess_li_pct_max = 0.05

lowaccess_pct_min = 0
lowaccess_pct_max = 0.05

priority_min = 0
priority_max = 10

def style_function_snap_rate(feature):
    return {
        'fillColor': colormap_snap_rate(feature['properties']['snap_rate']),
        'color': 'black',
        'weight': 1,
        'fillOpacity': 0.7
    }

def style_function_chd_pct(feature):
    return {
        'fillColor': colormap_chd_pct(feature['properties']['chd_pct']),
        'color': 'black',
        'weight': 1,
        'fillOpacity': 0.7
    }

def style_function_lowaccess_li_pct(feature):
    return {
        'fillColor': colormap_chd_pct(feature['properties']['lowaccess_li_pct']),
        'color': 'black',
        'weight': 1,
        'fillOpacity': 0.7
    }

def style_function_lowaccess_pct(feature):
    return {
        'fillColor': colormap_chd_pct(feature['properties']['lowaccess_pct']),
        'color': 'black',
        'weight': 1,
        'fillOpacity': 0.7
    }

def style_function_spa(feature):
    return {
        'color': spa_colors.get(feature['properties']['SPA_NAME'], 'white'),  # Boundary color
        'weight': 1,
        'fillOpacity': 0,
        "dashArray": "5, 5"
    }

def style_function_priority(feature):
    return {
        'fillColor': index_colors.get(feature['properties']['priority_index'], 'white'),
        'weight': 0,
        'fillOpacity': 0.7
    }

colormap_snap_rate = cm.linear.YlOrRd_09.scale(snap_rate_min, snap_rate_max).to_step(n=5)
colormap_chd_pct = cm.linear.Blues_09.scale(chd_pct_min, chd_pct_max).to_step(n=5)
colormap_lowaccess_li_pct =cm.linear.YlOrRd_09.scale(lowaccess_li_pct_min, lowaccess_li_pct_max).to_step(n=5)
colormap_lowaccess_pct = cm.linear.Blues_09.scale(lowaccess_pct_min, lowaccess_pct_max).to_step(n=5)

colormap_snap_rate.caption = "SNAP Rate"
colormap_chd_pct.caption = "CHD Percentage"
colormap_lowaccess_li_pct.caption = "Low Access Low Income Percentage"
colormap_lowaccess_pct.caption = "Low Access Percentage"

bounds = gdf.total_bounds

center_lat = (bounds[1] + bounds[3]) / 2
center_lon = (bounds[0] + bounds[2]) / 2

m = folium.Map(location=[center_lat, center_lon], zoom_start=8)

geojson_data = gdf.to_json()

#snap_rate_layer = folium.GeoJson(
#    geojson_data,
#    style_function=style_function_snap_rate,
#    name='SNAP Rate',
#    popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'], 
#                              aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'], 
#                              labels=True)
#).add_to(m)

chd_pct_layer = folium.GeoJson(
    geojson_data,
    style_function=style_function_chd_pct,
    name='CHD Percentage',
    popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'], 
                              aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'], 
                              labels=True)
).add_to(m)

#lowaccess_li_pct_layer = folium.GeoJson(
#    geojson_data,
#    style_function=style_function_lowaccess_li_pct,
#    name='Low Income Low Access Percentage',
#    popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'], 
#                              aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'], 
#                              labels=True)
#).add_to(m)

#lowaccess_pct_layer = folium.GeoJson(
#    geojson_data,
#    style_function=style_function_lowaccess_pct,
#    name='Low Access Percentage',
#    popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'], 
#                              aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'], 
#                              labels=True)
#).add_to(m)

spa_layer = folium.GeoJson(
    geojson_data,
    style_function=style_function_spa,
    name='SPA Overlay',
    popup=folium.GeoJsonPopup(fields=['SPA_NAME'], aliases=['Service Planning Area:'], labels=True)
).add_to(m)

index_layer = folium.GeoJson(
    geojson_data,
    style_function=style_function_priority,
    name='FI Index Overlay',
    popup=folium.GeoJsonPopup(fields=['priority_index'], aliases=['Priority Index:'], labels=True)
).add_to(m)

#colormap_snap_rate.add_to(m)
colormap_chd_pct.add_to(m)
#colormap_lowaccess_li_pct.add_to(m)
#colormap_lowaccess_pct.add_to(m)

folium.LayerControl().add_to(m)

m
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook